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ABSTRACT 

Quantum  mechanical  calculations  based  on  density 
functional  theory  (DFT)  are  used  to  study  dynamic 
behavior  of  shocked  energetic  materials  (EM).  In  this 
work,  we  present  results  of  quantum  molecular  dynamics 
(QMD)  simulations  of  shocked  pentaerythritol  tetranitrate 
(PETN),  a  conventional  high  explosive,  and  the  polymeric 
cubic  gauche  phase  of  nitrogen  (cg-N),  proposed  as  an 
environmentally  acceptable  energetic  alternative  to 
conventional  explosive  formulations.  These  simulations, 
made  possible  through  a  Challenge  grant  awarded  by  the 
DoD  High  Performance  Computing  Modernization 
Program  (DoD  HPCMP),  represent  the  leading  edge  of 
DFT  simulations  in  both  system  size  and  simulation  time 
with  over  4000  atoms  and  up  to  ten  thousand  time  steps 
utilizing  as  many  as  512  processors  per  run. 

1.  INTRODUCTION 

The  changing  nature  of  international  conflicts 
requires  the  transformation  of  war-fighting  capabilities 
within  the  Department  of  Defense.  The  attainment  of  the 
transformation  of  the  military  force  structure  depends  on 
the  development  of  lighter-weight,  more  robust  and 
integrated  armor  materials  and  for  systems  to  meet  the 
requirements  for  speed,  mobility  and  rapid  deployment  of 
weapon  platforms  and  to  enhance  combat  personnel 
effectiveness  and  survivability  in  the  face  of  a  widening 
range  of  more  lethal  and  countering  threats.  New 
energetic  materials,  with  substantially  enhanced 
performance  and  reduced  vulnerabilities,  must  be 
developed  for  future  weapons  systems  to  meet  these 
criteria.  Advanced  Energetic  Materials  (AEM)  are 
required  to  enable  high  priority  military  missions  ranging 
from  Hard  and  Deeply  Buried  Target  Defeat,  to  Advanced 
Propulsion,  to  lightened  highly  mobile  force  evolution 
and  the  thrust  towards  miniaturized  munitions  and 
systems.  It  is  recognized  that  weapons  superiority  is 
dependent  on  the  development  of  AEM  since  they  are  key 
to  enhanced  lethality  and  survivability  in  logistical  and 
tactical  environments.  Current  developmental  strategies 
for  formulation  of  advanced  EMs  focus  on  either  (1) 
tailoring  the  dynamic  response  of  existing  EMs  that 
would  permit  optimal  performance  with  minimal 
vulnerability  for  a  wide  range  of  future  gun,  missile, 
warhead,  active  protection  systems  and  emerging 
weapons  concepts;  or  (2)  increasing  energy  content  in 


candidate  materials  by  introducing  or  increasing 
molecular  strain  within  the  molecules.  With  regard  to  the 
first  strategy,  it  is  imperative  that  the  fundamental 
mechanisms  controlling  conversion  of  the  material  to 
product  be  completely  understood;  such  information 
could  then  be  exploited  to  tailor  the  material  or  weapon  to 
control  the  dynamic  response.  In  addressing  the  second 
strategy,  in  which  energy  content  is  increased  by 
introducing  strain,  conventional  strategies  currently 
applied  to  existing  EMs  must  be  discarded  in  favor  of 
completely  new  types  of  materials  due  to  limitations  on 
the  amount  of  strain  that  can  be  imposed  without 
destabilizing  the  system  [Dlott,  2004].  One  potential  new 
class  of  EM  involves  high  pressure  phases  of  covalent 
solids,  in  which  structural  energy  is  trapped  in  metastable 
states.  If  the  stored  structural  potential  energy  can  be 
liberated  quickly  enough,  it  is  possible  that  explosion  can 
occur  with  energies  several  orders  of  magnitude  larger 
than  conventional  explosives.  Subsequently,  these 
materials  could  provide  enhanced  or  configurable 
performance  in  advanced  weapons  applications. 

To  explore  these  emerging  candidate  materials  and  to 
obtain  a  fundamental  understanding  of  the  response  of  an 
EM  to  shock,  we  have  performed  quantum  mechanical 
calculations  based  on  density  functional  theory  (DFT)  to 
study  dynamic  behavior  of  both  conventional  and  novel 
shocked  energetic  materials  (EM).  In  this  work,  we 
present  results  of  quantum  molecular  dynamics  (QMD) 
simulations  of  shocked  pentaerythritol  tetranitrate 
(PETN),  a  conventional  high  explosive,  and  shocked 
polymeric  nitrogen  (cubic  gauche,  cg-N)  proposed  as  an 
environmental  acceptable  energetic  alternative  to 
conventional  EMs.  While  nitrogen  seems  far-removed 
from  conventional  EMs  (which  are  typically  large, 
polyatomic  organic  molecular  crystals,  as  seen  in  our 
depiction  of  the  simulation  cell  used  for  PETN  in  Fig.  1), 


Figure  1.  Snapshot  of  the  equilibrated  PETN  system 
used  in  the  QMD  simulations. 
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under  pressure  it  has  been  shown  to  take  on  a  polymeric 
form.  The  cg-N  form  is  a  bulk  three  dimensional  covalent 
network  like  diamond,  in  which  every  atom  has  three 
neighbors  compared  to  diatomic  molecular  forms  of 
nitrogen  (Fig.  2). 


Figure  2.  Snapshot  of  the  equilibrated  cg-N  system  used 
in  the  QMD  simulations. 

Since  diatomic  nitrogen  has  the  strongest  bond  in 
nature,  the  amount  of  energy  released  from  the  polymeric 
form  as  it  transitions  to  the  more  stable  standard  state  is 
extremely  large.  The  products  of  the  conversion  are  very 
simple;  the  solid,  which  is  composed  of  singly-bonded 
nitrogen,  transitions  to  the  diatomic  triply-bonded  gaseous 
form  and  small  chains.  The  amount  of  energy  available 
for  release  through  complete  conversion  is  several  times 
that  of  the  conventional  explosive  cyclotrimethylene- 
trinitramine  (RDX)  by  weight  and  volume.  Simulation 
has  confirmed  that  cg-N  is  stable  at  low  temperatures  and 
pressures,  but  will  react  very  rapidly  and  release  large 
quantities  of  energy  when  sufficiently  perturbed  [Mattson, 
2003].  The  simplicity  and  reactivity  of  this  single-element 
system  make  it  an  ideal  physically  realistic  test  case  for 
simulation  of  shock  initiation  in  a  material  trapped  in  a 
high-pressure  metastable  state.  In  addition  to  presenting 
results  for  shocked  conventional  and  notional  EMs,  this 
paper  will  briefly  describe  the  theoretical  methods,  the 
computational  details,  software  explored  to  enable 
adequate  treatment  of  the  calculations,  and  results  of  the 
calculations  to  date. 


2.  COMPUTATIONAL  DETAILS 

DFT  is  not  considered  the  most  accurate  QM  method 
available,  but  it  has  become  a  state-of-practice  in  the 
scientific  community  due  to  its  computational  efficiency 
relative  to  other,  more  accurate  methods.  For  example, 
the  CCSD(T)  [Kimmel,  2002]  quantum  chemistry 
techniques  are  formally  0(N7)  where  N  reflects  the 
number  of  electrons  in  the  system  in  sharp  contrast  to 
DFT  methods  which  typically  range  from  O(N)  -  0(N3). 
In  the  DFT  formalism  each  electron  is  treated 
independently  and  interacts  with  other  electrons  not 
directly  but  only  through  an  effective  potential  which  is  a 
function  of  the  total  charge  density.  This  mean-field 
approximation  drastically  reduces  the  computational 
effort  compared  to  higher  level  quantum  chemistry 
methods  while  maintaining  accuracy  sufficient  for  many 
applications.  Advances  in  approximations,  efficient 
mathematical  algorithms,  and  implementations  on 


scalable  platforms  have  further  increased  application  of 
DFT  to  large  and  complex  systems  of  interest  to  the  DoD. 
In  this  study,  we  have  utilized  the  PBE  [Perdew  et  al., 
1996]  form  of  the  Generalized  Gradient  Approximation 
(GGA)  of  the  DFT  as  implemented  in  the  local  orbital 
basis  code  CP2K.  For  nitrogen  we  have  used  the  double- 
zeta  valence  polarization  basis  (DZVP)  which 
corresponds  to  a  total  of  thirteen  basis  functions  per 
nitrogen  atom  Total  energies  and  forces  were  converged 
to  1.6  x  10"4  and  2.6  x  10"5  a.u.,  respectively,  and  were 
used  throughout  these  calculations. 

Two  sets  of  simulations  were  performed.  The  first  is 
a  DFT  molecular  dynamics  (MD)  simulation  of  a  filament 
of  material  subjected  to  flyer-plate  impact.  This  is  a  non¬ 
equilibrium  molecular  dynamics  (NEMD)  simulation  in 
which  the  dynamic  response  of  the  material  can  be 
directly  observed  through  monitoring  the  time  progression 
of  atomic  positions  and  velocities.  These  simulations 
were  used  to  study  PETN  and  cg-N.  The  second  set  also 
utilizes  DFT-MD,  except  the  simulations  calculate 
properties  of  the  material  under  thermodynamic 
equilibrium  at  various  conditions.  Equation-of- State 
(EOS)  information  is  obtained  from  these  simulations  and 
used  to  generate  the  shock  Hugoniot,  a  special  set  of 
points  satisfying  the  conservation  equations  of  mass, 
momentum  and  energy  across  the  shock  front  between  the 
quiescent  crystal  and  its  final  shocked  state.  This 
information  will  be  generated  using  the  Erpenbeck 
[Erpenbeck,  1992]  approach,  named  after  an  author  who 
used  classical  MD  to  generate  the  shock  Hugoniot  of  a 
simple  model  material.  The  shock  Hugoniot  was 
calculated  only  for  cg-N,  due  to  inherent  deficiencies  in 
DFT  to  adequately  describe  conventional  molecular 
organic  crystals  (such  as  PETN)  at  low  compression 
[Byrd  et  al,  2004].  These  limitations  in  DFT  do  not  affect 
the  cg-N,  a  covalently  bound  system. 

2.1.  Direct  Shock  Simulations 

In  both  PETN  and  cg-N  simulations,  the  simulation 
cells  are  composed  of  the  material  in  an  equilibrated 
configuration.  The  simulation  cell  for  cg-N  consists  of  a 
filament  of  energetic  material  composed  of  2304  atoms 
arranged  in  the  equilibrium  cg-N  configuration  at  T=250 
K,  P  =  1  atm.  The  initial  simulation  cell  for  PETN  is 
composed  of  an  8  x  2  x  2  block  of  unit  cells,  with  the 
molecules  arranged  in  the  experimental  configuration  and 
atomic  velocities  assigned  to  correspond  to  room 
temperature.  Because  DFT  does  not  adequately  describe 
PETN  at  the  ambient  state,  we  equilibrated  the  simulation 
cell  at  a  fixed  volume  corresponding  to  the  experimental 
values  at  room  conditions.  Atomic  velocities  were 
initially  assigned  to  the  atoms  to  correspond  to  room 
temperature.  All  other  degrees  of  freedom  were  allowed 
to  relax.  For  the  shock  simulations,  periodic  boundary 
conditions  are  imposed  in  the  two  directions 
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perpendicular  to  the  direction  of  shock  propagation.  In 
this  study,  the  shock  wave  will  initiate  at  the  far  left  edge 
of  the  filament  (see  Figs.  1  and  2)  and  propagate  from  left 
to  right.  Initial  velocities  consistent  with  the  desired  flyer 
plate  impact  velocity  are  assigned  to  the  atoms  within  a 
small  segment  of  the  material  at  the  far  left  edge  of  the 
filament.  The  remaining  atoms  in  the  filament  are  at  a 
thermodynamic  equilibrium,  with  the  exception  of  a  small 
rigid  block  of  atoms  at  the  far-right-edge  of  the  simulation 
cell.  This  is  to  ensure  that  the  surface  effects  at  the  free 
edge  will  not  influence  the  shock  propagation.  The 
resulting  integration  of  the  microcanonical  equations  of 
motion  using  DFT-MD  for  this  set  of  initial  conditions 
will  simulate  shock  initiation  through  flyer-plate  impact. 
The  position  of  the  shock  discontinuity  through  the 
filament  is  monitored,  and  when  it  approaches  the  end  of 
the  filament  (i.e.  the  atoms  that  are  fixed),  additional 
material  will  be  inserted  before  the  fixed  atom  region, 
ahead  of  the  progressing  shock  wave.  The  state  of  the 
material  that  is  added  is  also  in  thermodynamic 
equilibrium  and  its  atomic  arrangements  are  consistent 
with  that  of  unshocked  material  at  that  thermodynamic 
state. 

Since  a  shock  wave  travels  faster  than  the  sound 
speed  of  the  material,  there  is  no  need  to  simulate  large 
portions  of  undisturbed  material  ahead  of  the  shock  front; 
rather,  we  have  found  in  earlier  classical  MD  studies  that 
adding  material  ahead  of  the  propagating  shock  front  is 
sufficient  to  adequately  describe  the  dynamics  of  the 
process.  This  method  has  been  used  in  many  classical 
MD  simulations  of  shock  waves  [Rice,  1998].  Obviously, 
in  such  a  simulation  scheme,  the  simulation  size  grows  as 
the  shock  wave  progresses  through  the  material, 
significantly  increasing  the  computational  demands 
during  the  simulation. 

2.2  Evaluation  of  the  Shock  Hugoniot 

The  thermodynamic  quantities  of  a  material  in  the 
quiescent  and  final  shocked  states  are  related  by  the 
conservation  equations  of  mass,  momentum,  and  energy 
across  the  shock  front,  and  are  represented  by  the  shock 
Hugoniot,  a  set  of  EOS  points  that  satisfy  the  Hugoniot 
function: 

Hg(r,  v)  =  0  =  E-Eo- y2(P  +  p0)(v0  -V)  (i) 

In  this  equation,  E,  P,  and  V  are  the  internal  energy/unit 
mass,  pressure,  and  volume/unit  mass,  respectively.  E0, 
P0,  and  Vo  are  the  internal  energy/unit  mass,  pressure  and 
volume/unit  mass  in  the  quiescent  state.  The  procedure 
for  identifying  points  on  the  shock  Hugoniot  using  the 
method  of  molecular  dynamics  is  straightforward 
[Erpenbeck  1992].  For  each  volume,  a  series  of  DFT-MD 
simulations  of  the  NVT  ensemble  are  performed  over  a 
range  of  temperatures.  Once  convergence  of  the 


thermodynamic  properties  is  reached  for  each  simulation, 
the  Hugoniot  function  is  evaluated.  The  resulting  set  of 
calculations  will  produce  a  series  of  Hg(7)F,  which  is  then 
fitted  to  a  polynomial  from  which  the  Hugoniot 
temperature  (rHg)  can  be  identified  from  Hg(rHgV  =0. 
Once  this  has  been  done,  the  dependence  of  either  the 
pressure  or  the  internal  energy  on  temperature  can  be 
determined  through  fitting  these  quantities  to  polynomials 
and  evaluating  them  at  T\\g-  Each  of  the  calculations  used 
to  construct  the  shock  Hugoniot  for  cg-N  was  composed 
of  256  atoms,  with  their  initial  arrangement  consistent 
with  that  of  the  cg-N  structure  at  T=298K,  P  =  1  atm. 
Each  trajectory  was  first  equilibrated  for  4000  time  steps 
and  thermodynamic  values  were  averaged  over  the 
subsequent  4000  time  steps  (1  time  step=  1  fs). 


3.  RESULTS 

Two  types  of  MD  calculations  were  performed  to 
provide  insight  into  the  nature  of  shocked  cg-N.  The  less 
computationally  demanding  is  the  generation  of  the  shock 
Hugoniot  using  equilibrium  MD  methods.  This 
methodology  does  not  require  large  numbers  of  atoms, 
and  uses  MD  to  generate  the  EOS  for  the  materials.  Its 
computational  cost  comes  from  the  number  of  EOS  points 
needed  to  map  out  the  shock  Hugoniot.  Each  point  on  the 
shock  Hugoniot  was  obtained  by  interpolating  the 
Hugoniot  function  for  four  EOS  points.  In  total  this  has 
required  in  excess  of  500,000  CPU  hours  on  the  SGI 
Origin  3900  at  AFRL  MSRC. 

As  evident  in  Fig.  3,  no  anomalous  behavior  was 
found  in  the  principal  shock  Hugoniot  for  cg-N. 
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Figure  3.  Principal  shock  Hugoniot  of  cg-N  in  pressure- 
density  (left)  and  pressure-temperature  (right)  planes. 

However,  additional  NVT  simulations  showed  that  cg-N 
will  undergo  a  phase  transition  at  approximately  P  =  50 
GPa,  T  =  4000K;  the  material  transitioned  from  the 
polymeric  solid  form  to  liquid  composed  of  both  atoms 
and  molecules.  The  material  in  this  state  is  not  well 
characterized.  Figure  4  shows  three  snapshots  of  the 
system  at  various  time  steps  during  the  trajectory 
integration.  The  left-most  frame  shows  the  system  in  the 
initial  cg-N  arrangement,  and  the  right-most  frame  shows 
a  distribution  of  atoms  and  clusters  ranging  in  size  from 
two  to  ten  atoms.  The  structure  of  the  material  is 
consistent  with  a  liquid. 
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Figure  4.  Snapshots  of  cg-N  undergoing  a  phase 
transition  at  t=0  (left-most  frame),  1=0.35  ps  (middle 
frame)  and  1=0.7  ps  (right-most  frame)  at  T=4000  K, 
P=50  GPa. 

The  dynamic  behavior  of  the  material  during  this  portion 
of  the  trajectory  shows  that  the  clusters  are  not  stable  as 
the  atoms  rapidly  change  partners.  It  is  not  known 
whether  any  true  cluster  formation  is  obtained.  The 
middle  frame  is  representative  of  a  transient  state  of  the 
material  that  appears  to  have  some  structure.  Further 
calculations  are  needed  in  order  to  make  a  definitive 
characterization  of  the  transient  and  final  states  of  the 
material. 

Our  first  shock  simulations  were  performed  on  the 
polymeric  form  of  nitrogen.  We  have  simulated  shock 
wave  propagation  through  the  cg-N.  The  mechanical 
shocks  were  initiated  by  impacting  the  crystal  with  a  low 
flyer  plate  initial  velocity  to  determine  the  influence  of  a 
mild  perturbation  to  this  metastable  state.  Simulations  in 
which  periodic  boundary  conditions  are  imposed  in  all 
dimensions  were  performed;  however,  in  one  type  of 
simulation  a  vacuum  is  placed  behind  the  flyer  plate.  In 
the  other  type,  no  vacuum  is  present.  In  both  types  of 
simulations  a  small  segment  of  crystal  (54  atoms)  on  the 
opposite  end  of  the  crystal  are  held  fixed.  The  results  of 
shock  wave  simulations  using  the  two  types  of  simulation 
cells  (i.e.  with  or  without  a  vacuum  behind  the  flyer  plate) 
produced  different  results.  For  periodic  systems  without  a 
vacuum  behind  the  flyer  plate  or  defects,  low-velocity 
impacts  were  insufficient  to  produce  reaction  within  the 
crystal.  Rather,  the  shock  wave  merely  traversed  the 
filament  between  the  flyer  plate  and  the  rigid  material  at 
the  opposite  side  of  the  cell  several  times,  with  no 
reactive  events.  However,  reactive  events  were  observed 
in  the  simulations  in  which  a  vacuum  was  present  behind 
the  flyer  plate  (Fig.  5). 

However,  this  reaction  does  not  occur  upon  shock  wave 
compression  and  subsequent  expansion.  After  reaching 
the  rigid  slab  of  molecule,  the  shock  was  reflected  back 
through  the  shocked  material.  Upon  reaching  the  left¬ 
most  edge  of  the  simulation  cell,  the  material  at  the  free 
edge  of  the  crystal  unraveled  in  chains  of  up  to  ten 
nitrogen  atoms.  These  chains  subsequently  decompose 
into  azide  and  diatomic  nitrogen  molecules,  as  shown  in 
Fig.  5. 
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Figure  5.  Snapshot  of  the  free  edge  of  the  shocked  cg-N 
at  t=3.6  ps. 

QMD  simulations  of  the  filament  with  the  free  edge 
under  equilibrium  conditions  (i.e.  no  shock,  250  K,  1  atm) 
does  not  produce  unraveling  of  the  free  edge,  although  the 
filament  undergoes  an  expansion  and  the  atoms  at  the 
edge  are  more  mobile.  It  appears  that  the  unraveling  of 
the  material  and  its  subsequent  decomposition  into  tri- 
and  diatomic  molecular  species  is  a  result  of  the  reflected 
shock  wave.  We  are  continuing  further  simulations  in 
which  the  shock  wave  is  not  allowed  to  reflect. 

We  are  repeating  the  same  procedure  for  the 
conventional  energetic  material  pentaerythritol  tetranitrate 
(PETN).  We  have  completed  the  equilibration,  and  are 
running  the  shock  wave  simulations.  The  four  snapshots 
in  Fig.  6  show  the  system  at  the  beginning  of  the 
trajectory  (top-most  frame),  0.109  ps  (second  frame),  0.28 
ps  (third  frame),  and  0.418  ps  (bottom  frame)  during  the 
trajectory. 

To  date,  the  shock  wave  has  traversed  the  filament 
and  reached  the  rigid  slab  of  molecules  at  the  far-right 
edge.  PETN  is  considerably  more  complex  than  cg-N  due 
to  the  number  of  degrees  of  freedom  and  types  of  reaction 
products.  The  structural  nature  of  PETN,  which  is  a 
weakly  bound  molecular  crystal  in  which  energy  release 
is  due  to  chemical  reaction,  differentiates  it  from  cg-N,  a 
covalent  solid  in  which  the  energy  release  is  due  to  a 
phase  transition.  Therefore  it  is  not  surprising  to  see  that 
the  reaction  mechanisms  described  for  cg-N  are  not  the 
same.  These  differences  are  apparent  in  the  series  of 
snapshots  shown  in  Fig.  6.  In  the  second  frame  (t=0.109 
ps),  the  shock  wave  has  started  to  traverse  the  material, 
and  the  material  behind  the  shock  front  has  begun  to 
expand.  In  the  third  frame  (t=0.28  ps),  the  region  through 
which  the  shock  wave  has  passed  appears  compressed  and 
disordered,  and  reaction  products  such  as  NO2  and  OH,  as 
well  as  some  radical  fragments  of  the  PETN  skeleton,  are 
evident.  Finally,  at  t=0.418  ps,  we  see  a  material  profile 
that  suggests  rarefaction  in  the  density;  reactions  are 
clearly  evident  in  this  region. 
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represent  the  largest  QMD  simulations  on  explosives  to 
date. 


Figure  6.  Snapshots  of  the  PETN  system  at  t=0  (top¬ 
most  frame),  0.109  ps  (second  frame),  0.28  ps  (third 
frame),  and  0.418  ps  (bottom  frame). 

To  our  knowledge,  these  are  the  largest  DFT  simulations 
of  an  energetic  material  that  have  ever  been  performed. 


CONCLUSIONS 

Quantum  MD  simulations  have  been  used  to  provide 
an  atomic-level  description  of  the  material  response  of  a 
conventional  explosive,  PETN  and  a  novel  candidate 
energetic  material,  a  polymeric  form  of  solid  nitrogen 
known  as  cg-N,  to  shock.  Both  simulations  showed 
reactions  after  passage  of  the  shock  wave;  however,  the 
mechanistic  details  differed.  PETN  reaction  occurs 
through  dissociation  of  N02;,  OH,  and  transient  reactive 
molecular  radicals,  whereas  the  cg-N  material  at  the  free 
edge  of  the  shocked  material  unraveled  in  chains  that 
subsequently  decomposed  into  di-  and  tri-nitrogen 
molecules.  Additionally,  we  calculated  the  shock 
Hugoniot  of  cg-N,  for  future  comparison  with  results 
from  direct  shock  simulations.  These  calculations 
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